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Abstract 


This paper aims to apply and investigate the compact finite difference 
methods for solving integer-order and fractional-order Riccati differential 
equations. The fractional derivative in the fractional case is described in 
the Caputo sense. In solving the Riccati equation, we first approximate 
first-order derivatives using the approach of compact finite difference. In 
this way, the system of nonlinear equations is obtained, which solves the 
Riccati equation. In addition, we examine the convergence analysis of the 
proposed approach for the fractional and nonfractional cases and prove that 
the methods are convergent under some suitable conditions. Examples are 
also given to illustrate the efficiency of our method compared to other 
methods. 
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1 Introduction 


In recent years, there has been a growing interest in fractional computation 
[14, 26, 31, 33, 34]. Fractional differential equations have become increasingly 
important as they have applications in various fields of science and engineer- 
ing [13]. Numerous phenomena in fluid mechanics, viscoelasticity, chemistry, 
physics, finance, and other sciences can be described successfully by mod- 
els using mathematical tools of fractional calculation, that is, the theory of 
fractional-order derivatives and integrals. Much important work on theoreti- 
cal analysis [38, 10] has been carried out, but the analytical solutions of most 
fractional differential equations cannot be achieved explicitly. Numerical so- 
lution strategies based on convergence and stability analysis were used by 
many authors [12, 11, 13, 16, 20, 35, 36, 39, 41, 22]. Liu has carried out 
extensive research on the finite difference method of fractional differential 
equations [22, 23, 24]. The two most frequently used are the Riemann-— 
Liouville and Caputo type. The difference between the two definitions is in 
the order of evaluation [29]. 


In this paper, we consider the following Riccati equation: 


u(x) = p(x) + q(x)u(x) + r(x)u? (2), 0<a<T, (1) 
u(0) = 0. 


Also, we consider the following fractional Riccati equation: 


D%u(x) = p(x) + ¢(x)u(x) + r(x)u?(2), 0<a<l, 0<a<T, (2) 
along with the initial condition 
u(0) = 0, (3) 


where x € R and p(x), q(x), and r(x) are known functions. Moreover, D® is 
the Caputo derivative operator of the fractional-order a, which is defined as 
below: 


1 xv 
D°u(x) = Tad —a) | (x — 8) °u'(s) ds. (4) 
In the past, two scholars, Bernoulli (1654-1705) and Riccati (1676-1754) in- 
troduced and assessed a particular case of differential equations (2). The 
Riccati differential equations (RDEs) and fractional Riccati differential equa- 
tions (FRDEs) are used in many physical phenomena. Such applications can 
include control systems, robust stabilization, diffusion problems, network 
synthesis, optimal filtering, stochastic theory, controls, financial mathemat- 
ics, optimal control, river flows, robust stabilization, financial mathematics 
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dynamic games, linear systems with Markovian jumps, stochastic control, 
econometric models, and invariant embedding [32, 28, 4, 19, 15, 9, 21, 5, 30]. 
Many researchers have used numerical approaches to solve the RDEs and 
FRDEs. Some standard procedures can be referenced, including the differen- 
tial transform method [7], series solutions Adomian’s decomposition method 
[1], Homotopy perturbation method [1], variational iteration method [18], 
Homotopy analysis method [37], piecewise spectral-collocation method [6], 
and so on [8, 27, 25, 3]. 

This paper aims to obtain numerical solutions to (1)—(3) using a high-order 
compact finite difference approach. 


Several researchers have employed the compact finite difference method 
to solve fractional differential equations. Du, Cao, and Sun [14] have used 
the compact finite difference method to solve the fractional diffusion-wave 
equation. Gao and Sun [17] have also employed the compact finite differ- 
ence method to solve the fractional sub-diffusion equation. They have also 
proved the stability and convergence of their method. Cui [13] solved the 
one-dimensional fractional diffusion equation via a high-order compact finite 
difference scheme and obtained a fully discrete implicit system by Grunwald— 
Letnikov’s discretization of the Riemann—Liouville derivative. 


The present study is organized as follows: In Sections 2, 3, and 4, the 
compact finite difference methods are reviewed and applied to solve (1)- 
(3). Also, their convergence is discussed. In Section 5, the numerical results 
obtained by the proposed methods are presented. We also compare the results 
of our approach and those of the proposed methods in [2]. The conclusion 
and the advantages of the proposed technique are presented in Section 6. 


2 Compact finite difference scheme 


In this work, our primary goal is to apply the compact finite difference method 
to solve (1)—(3). For this, we first subdivide the range 0 < « < T to N equal 
partitions with step length h as follows: 
=0 ; = th i= 0,1 N a (5) 
ty =V, Ly = th, ’=U,1,...,4V, = 
Set 


uj © u(x), ti. PSU (Bi) 


For the first derivatives, the following compact finite difference scheme was 
given in [40]: 
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/ {x= pe 4 1 
at + Us, = n (79 uo : Adu, + 6u2 3 U3 1b U4), 
uy_y + 4uy + Uj, = F(—t-1 + Ui41), i=1,...,N-—1, (6) 
wy 9 t4uy_, = 2(—4u + 4u — 6u + 4u + tun) 

N-2 N-1”7 h 4 N-4 3 N-3 N-2 N-1 12 N)>+ 


All above relations have the accuracy of O(h*). The matrix form for (23) 


is 
posi 
Aju = —Byu, (7) 
h 
where 
04 1 0...0 
14 1 0...0 
A, =| 0 - 0 ; 
bacco: At td al 
er 0 ; i (N-+1)x (N41) 
~ 727 6 3. a 0 ... O 
-3 0 3 0 00... 0 
0-30 3 00... 0 
BE) 3G 0 
>: ... O 0 0-3 0 3 
1 4 11 
Ob ees Oa a Oh aoe (N41) x(N-41) 
Also, u = [uo, u1,---, un]? and u’ = [up, ul,...,uy]?. 


Lemma 1. The coefficient matrix A, is invertible. 


Proof. Let us expand A, along the first column. Then 


41 0...0 

14 1. 
det(A,) = —det | : ° a4 

Q...1 41 

Ose ds A Oe eng 


det(A,) =(—-1)%det | i+. +... : 40. 
4 1 
14 (N—1)x(N—1) 
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According to Lemma 1, from (24), we have u’ = ZA, Biu. By defining 
C = A;'By, the following relation holds for wu’: 


ia > Cu, (8) 


and in the component form, we have 
1x 
uj; = pests ¢=0,...,N. (9) 


Lemma 2. The coefficient matrix B, is invertible. 


Proof. Let us expand B, along the first row. Then 


det(B,) =—det | 0 
0 


2 ate 0 
0... 0 - 


1 
4 


0 
...0 

det(By)=(-1)Ndet | 9 ee tg £0. 
3 


SW OD 8 O58 8 singe, 


According to Lemma 2, it follows that the matrix C is invertible. 


3 Compact finite difference scheme for Riccati problem 
in a@ = 1 case and its convergence 


This section uses the compact finite difference scheme for the nonfractional 
Riccati problem and investigates its convergence. Consider the subsequent 
classical Riccati initial value problem 


u(x) = p(x) + ¢(x)u(x) + r(x)u? (x), 0<a<T. (10) 
Its initial condition is 


u(0) = 0. (11) 
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Using (10), we have 
Up = (xo). (12) 


So, using (9), equation (12) can be written as 


N 
me C1,3+1Uj = P(Xo). (13) 


A 


For « = 2;, one can write (10) as 
u' (a4) = p(ai) + g(a) ula) + r(ai)u? (2), (HM IN. (14) 


Thus from (9) 


Re 


N 
i oy Cag ade — Dey) — gas us = r(a,)u? = 0, t#=1,...,N. (15) 


Equations (13) and (15) form a system including N + 1 equations and N +1 
unknowns uo, Ui,.--,Un, that can be solved by Maple software. 


Now, the convergence analysis of the proposed method for (10) along with 
initial conditions (11) is investigated. 


Theorem 1. Let U = [u(xo), u(x1),-..,u(xn)]” be the vector of exact so- 
lution to (1) along with its initial condition, and let u = [uo,u,...,un]? be 
the numerical solution at the same points obtained by (13) and (15). Then 

||Ell < O(h7), (16) 


provided h||C~+||||M|| <1, where E = [eo,e1,...,ew]? and e; = u(a;) — us, 
i=0,...,N (||- || is the infinity norm). 


Proof. According to (13) and (15), for a numerical solution, we have 


N 
to ja0 C1.s+1Uj = P(2o), 


‘N ; 
eo jao Cit1.3+1Uy — P(e) — (asus — r(wa)uz = 0, t=1,...,N, 
and for an exact solution, we have 


{ tN ers 41u(ay) = p(xo) + O(h4), 
A peo C419 414(@y) — p(w) — a(wi)u(as) — r(wi)u? (a) = O(h4), 1=1,...,N. 


By subtracting (17) and (18), one concludes that 
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t pho €1,341(u(aj) — uy) = O(h*), 
Se ec tue 3) — Uj) — a(ai) (u(x) — ui) 
—r(a;)(u?(a;) — u?) = O(h4), i=1,...,N. 


I 


Using the Taylor expansion, we have 


2 
Oa) a? = Yaa (ule) — tu) FON), FTN. (20) 


In relation (19), we have 


i doje C1j41(u(aj) — uj) = O(h4), 
i ya ci41,j+1(u(xz) — uz) — a(xi)(u(ai) — ui) 
—2r(a;)u(a;)(u(a;) — us) = O(h*) + O(h?), $= 15 4053.N: 
(21) 


Thus, one concludes that 


otigere; = O(P), 
ear Ci+1,j +167 — ha(xi)ei — 2hr (ax; )u(a;)e; => O(h3), 1= Ly cat LN, 


(22) 
where e; = u(x;) — uj, 7 =0,...,N, and u; © u(a;). Therefore, (22) can be 
written as 

c11e0 + cize1 + c13€2 4 ci,n41en = O(h*), 
c21e90 + co2e1 + c23e2 +++: +c2,n41eEN — ha(ai)e1 — 2hr(a1)u(x1)e1 = O(h3), 
cg1eo + c32€1 + c33€2 4 c3,n+1eN — hq(w2)e2 — 2hr(a2)u(x2)e2 = O(h*), 
CN+1,1€0 + CN41,2€1 + CN 41,3€2 +°+*> +¢N41,N+41EN 
—hq(an)en — 2hr(2n)u(an)en = O(n). 
(23) 
The matrix form of the above equations is as follows: 
[(C—hQ-—hRJJE =T, (24) 
where Q = diag(0,q(a1),..-,¢(an)), R = diag(0,r(a1),...,r(aw)), J = 
diag(0, 2u(21) 2u(an)), and 
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O(n? 
O(h?) C11 C12 see C1,N41 
O he C21 C22 see) CQ2.N41 
T = | OCA") C= 
O(n?) EN+1 ENF12 ++ CNSLNTLS (4) x (N41) 


(N+1)x1 
(25) 


By replacing M = Q+ RJ in relation (24), we have [C —hAM]E = T. 
Because C’ is invertible, we can write 


(I-hC 'M)E=C"T. (26) 


Now, if we assume h||C~+||||M|| < 1, then we conclude the matrix I — 
hCc— iM is invertible. By the geometric series theorem, we have 


1 
| < : 
1—Al|C-3||||M|| 


(7 —nhCc-*M)-+ (27) 


From (26), we have EF = (I —hC7!M)"!C7!T. Thus |/E|] < || - 
hC~*M)“*"]“* IIIT 1- 


Now from relation (27), we can write ||E|| < rye OMIT. 


Because ||T'|| = O(h?), we can derive || E|| < — = O(h?). 


4 Implement the compact finite difference scheme for 
the fractional Riccati problem and its convergence 


In this section, we introduce a compact finite difference scheme for the frac- 
tional Riccati problem of order 0 < a < 1. According to (2), we rewrite the 


Caputo derivative in x = 2;,i=1,...,N,as 
Qos a (28) 
D*u(x;) = ——~ | ———ds. 28 
Lissa). (ae a)* 


Now, the above equation can be written as 
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k=0% 
i-1 
1 ; a = 
= ———u, (a; — s) “ds (29) 
T(1—-a) » Be 
— 1 ul. sti = t,)'% = (x; oe Lk41) = 
T(l-a) ' cam l-a 


T(l-a l-a 
hae (30) 
! a-1 
U, 
~ Fa-IP(2—q) Zw Wi 
ho-!T(2 — a) ame 


where aj;_, = (i— k)'~* — (ik —1)1-9, i= 1,...,N and kk =0,...,i-1. 


Thus, the solution to (2) can be approximated using the following equa- 
tions: 


i-l 
ui 2 . 
he-1T(2— a) 24 %-+ = p(ai)+q(a;)ur+r(aj)uZ, 0<a< 1, t= 1, sony N, 
(31) 
N 
where wu), = i YS ci41j41Uj, 1 = 1,...,N. In the matrix form, (31) is equiv- 
j=0 
alent to 
Ful = p(G+Qu+ Ru’), (32) 
where p = h°-!T(2—a), Q = diag(q(21),-.-,4(2N)), R= diag(r(21),...,7(en)), 
p(r1) 
p(x2) 
u=([u,...,un]?, Taal rare Plea G= . , and 
p(n) 
a, 0 0 0 0 
0 a, + ag 0 0 0 
0 0 a, +a2+ 43 0 0 
F=]| , ; ; : 
0 0 0 a, t+ag+---+ayn_1 0 
0 0 0 0 a, +a2+°::+an 


(33) 
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For i = 1,..., N, (31) can be used to form a system including N equations 
and N unknowns u1,...,uy, that can be solved by Maple software. 


Now, we discuss the issue of convergence. For convergence analysis of the 
fractional case, we need the following Lemma. 


Lemma 3. [35] Suppose u € C?(0, xi]. Then 


Law ee fo ni — 8) “ds| 


1 - ee Me 
12 2-—a 


(34) 


<i (1+2>9)] max ful"(s)|hPn*. 


From (30), we have 


i-1 


1 
S > a;_nu' (ai) + Ri, i= Ly-eba gl, (35) 


D%u(x;) = ———-~ 
ules) = Fe-IT@— a) var 


where according to Lemma 3 


_ 2—a 
R< 1 fe a 2 


ee (1+2-%)] max |u(s)|h?-*. (36) 


12 2-—a 0<s<x; 


For « = aj, by replacing (35) into (2), we have 


y a;—xu' (x2) = p(p(as) + a(as)Ju(as) +r@i)w(a)) +R, 1=1,...,N, 
k=0 


7 (37) 
where p = h°~'T(2—a) and R; = h®°'T(2—a)R;,i =1,...,N. 
In the matrix form, (37) is equivalent to 
FU’ = p(G+ QU + RU’) +R, (38) 
where F' is the matrix defined in relation (33), U’ = [u’(#1),...,u(aw)]’, 
U = [u(a1),...,u(zn)|*, and R= h?—11(2— @)[Ri,..., Rn]. 
Theorem 2. Let U = [u(z1),...,u(xy)|? be the vector of exact solution 
to (2) along with its initial condition at points x0,271,...,a., and let u = 
[ui,--.,Un]? be the numerical solution obtained by (31). Then 
||El| < O(h?-*), (39) 


provided h||C~+|||M + N|| <1, where EF =U —uand 
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2u(ai) 0 0 
0 2u(a2) 0 
J= ; 
0 ... 0 2u(ay) 


Proof. According to (38) and (32), for the exact and numerical solutions, we 
have 


FU’ = p(G+QU + RU?) +R, (40) 
Fu’ = p(G+ Qu+ Ru’). 
By using (40), one concludes that 
F(U’ — wu) = p(Q(U —u) + RU? — u?)) +.B. (41) 


Therefore, by replacing u’ = +Cu from (8) and U’! = +CU +T, into (41), 
we have 


1 


ZC —u) — pF QU —u)— pF CRU? -w)=FOR+T, (42) 


where T; = O(h*) is the local truncation error of system (23). 
Moreover, U? — u? can be written as 


u? (a1) — uy 
gwen | 4 | ean, (43) 
w2(ey) - uy 
where 
O(h?) 
O(h?) 
=| 
O(h?) 
Therefore, by replacing (43) into (42), we have 
~CE — pF ‘QE- pF 'RJE=F'R+T%+ pF RT. (44) 


By inserting relations M = pF~'Q and N = pF~'R into (44), it can be 
written as 


(C—hM —AN)E=h(F'R+T, + NT»), (45) 
(I—hC1(M+ N))E =hC1(F 1R4+T, + NT). (46) 
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Now, if Al|C7~|||| M+ N|| <1, then (I — hC~!(M + N) is invertible and 
E=h(I—hC-1(M + N))-!C71(F-'R+ T, + NT), 
|| Ell < |Z — hO*(M + NY) NOT ETA + Zl + IIZ2ID. 
It follows that 


h|}C*GLE- MRI + Zl) + VU Zall) 
1—Al|C~*|||[M + NI 


|E\| < , (47) 


_ Therefore, using the relations R = h°-!1P(2 — a)R and (36), we have 
|| R|| = O(A), so 
O(h?) | O(h?) 


|E|| < O(h?) + OG) 


= O(h?-*) 40(h®-%)+40(h3-*) = O(h?-°). 
(48) 


5 Numerical results 


This section applies our compact finite difference schemes to two examples to 
illustrate their effectiveness. Maple software is used for obtaining numerical 
results. 


Example 1. Consider the following fractional RDE as the first example: 


er 0<a<l, 0<a<T, 


u(0) = up = 0. ) 


exp(2x)—1 


The exact solution is u(a) = 

In Figure 1, a comparison between the exact solution for a = 1 and the 
numerical solution for a = 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1, and T = 1 is 
shown. Table 1 presents numerical solutions at some points of [0,1] and for 
different values of a, at T = 1. Table 2 presents a comparison between the 
exact solution for a = 1 and the numerical solution for T = 10. Also, Figure 
2 shows a comparison between the exact solution for a = 1 and the numerical 
solution for T = 10. 

We have calculated the rate of convergence of our methods (denoted by 
ROC) with the following formula: 


Error2? 


ROC = log,( i (50) 


Table 3 shows the obtained maximum errors and ROC for a = 1, T = 1, 
and N = 5,10, 20, 40, 80,160. Also, Figure 3 shows the numerical and exact 


Error’ 
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solutions for a = 1, T = 1, and N = 10. The numerical rate of convergence 
is highly consistent with our theoretical analysis results. 

In Table 4, we compare the approximate solution and exact solution of 
the present method with the trigonometric transform method (TTM) [2] at 
points 0.2, 0.4, 0.6,0.8,1, for a = 1. Also, in Table 5, we compare the error 
of solutions of the present method with TTM [2] for a = 1. 


Table 1: Exact solutions and numerical solutions of Example 1 for N = 10, 
T =1, and a = 0.8, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1 


a 0.3 0.6 0.7 0.8 0.9 0.95 0.99 0.999. 1 Exact 


0.1 | 5.38 x 10-1 | 2.66 x 10-! | 2.06 x 10-! | 1.60 x 10-1 | 1.25 x 10-1 | 1.11 x 107! J 1.01 x 10-1 | 9.98 x 10-? | 9.96 x 10-2 | 9.96 x 10-2 
0.2 | 7.43 x 10-1 | 4.34 x 10-1 | 3.53 x 10-1 | 2.88 x 107! | 2.37 x 10-1 | 2.16 x 10-1 | 2.00 x 107! | 1.97 x 10-1 | 1.97 x 10-1 | 1.97 x 107? 
0.3 | 8.34 x 10-1 | 5.51 x 10-! | 4.66 x 10-1 | 3.95 x 10-1 | 3.37 x 10-1 | 3.12 x 107! | 2.95 x 10-1 | 2.91 x 10-1 | 2.91 x 107! | 2.91 x 10-1 
0.4 | 8.83 x 10-1 | 6.38 x 10-1 | 5.57 x 10-1 | 4.86 x 10-1 | 4.27 x 107+ | 4.02 x 10-1 | 3.84 x 107! | 3.80 x 107! | 3.79 x 1071 | 3.79 x 107! 
0.5 | 9.14 x 10-1 | 7.05 x 10-1 | 6.31 x 107! | 5.64 x 1071 | 5.07 x 10-1 | 4.83 x 107! J 4.66 x 10-1 | 4.62 x 10-1 | 4.62 x 107! | 4.62 x 10-1 
0.6 | 9.34 x 10-! | 7.57 x 107! | 6.91 x 107! | 6.30 x 10-! | 5.79 x 10-! | 5.56 x 10-! | 5.40 x 107! | 5.37 x 10-! | 5.37 x 10-! | 5.37 x 10-1 
0.7 | 9.48 x 10-1 | 7.99 x 10-1 | 7.41 x 107! | 6.87 x 10-1 | 6.41 x 10-1 | 6.21 x 10-! J 6.07 x 10-1 | 6.04 x 10-1 | 6.04 x 107! | 6.04 x 10-1 
0.8 | 9.58 x 10-1 | 8.32 x 10-1 | 7.82 x 10-1 | 7.35 x 107! | 6.95 x 10-1 | 6.78 x 10-1 | 6.66 x 10-! | 6.64 x 10-1 | 6.64 x 10-1 | 6.64 x 107! 
0.9 | 9.65 x 10-1 | 8.59 x 10-! | 8.16 x 10-1 | 7.76 x 10-1 | 7.42 x 10-1 | 7.28 x 10-! J 7.18 x 10-1 | 7.16 x 10-1 | 7.16 x 107! | 7.16 x 10-1 
1.0 | 9.71 x 10-1 | 8.81 x 10-1 | 8.44 x 10+ | 8.10 x 10-1 | 7.82 x 10-1 | 7.70 x 10-1 | 7.63 x 107! | 7.61 x 10-! | 7.61 x 107! | 7.61 x 107-1 


Figure 1: Comparison between the exact solution of Example 1 for a = 1 
and numerical solutions for a = 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1 and T = 1 
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Table 2: Comparison between the exact solution and numerical solutions of 
Example 1 for a= 1, T = 10, and N = 100 


Numerical solution 


Exact solution 


Error 


0.7615917576 
0.9640223166 
0.9950446865 
0.9993096449 
0.9998709681 
0.9999133772 
0.9998538332 
0.9997188603 
0.9994538522 


OONDOBRWNYEH!S 


0.7615941559 
0.9640275800 
0.9950547536 
0.9993292997 
0.9999092042 
0.9999877116 
0.9999983369 
0.9999997749 
0.9999999695 


2.3983554 x 10-5 
5.2634336 x 107° 
1.0067173 x 107° 
1.9654751 x 1075 
3.8236123 x 107° 
7.4334413 x 1075 
1.4450373 x 1074 
2.8091453 x 1074 
5.4611733 x 1074 


Table 3: Maximum absolute errors and ROC of Example 1 fora = 1, T = 1, 


and N = 5, 10,20, 40, 80, 160 


N Maximum Absolute Error ROC 
5 7.95 x 10-4 - 
10 101 x10-? 5.38 
20 9.47 x 1077 4.33 
40 3.43 x 1078 4.79 
80 1.16 x 1079 4.88 
160 3.99 x 10714 4.87 
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1.2 
qe 
0.8 
06- 
0.4- 
0.27 —#—numerical solution} 7] 
= exact solution 
0 | | | | | | | | | 


0 1 2 3 4 5 6 fj 8 9 10 


Figure 2: Comparison between the exact solution and numerical solutions of 
Example 1 for a=1, T = 10, and N = 100 


1 T T 


—#*— Numerical solution 
- - -Exact solution 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
xX 


Figure 3: Comparison between the exact solution and numerical solutions of 
Example 1 for a=1, T=1, and N = 10 


Example 2. Let the following FRDE be the second example 
D%u(x2) = 1+ 2u(x) — u?(2), 0<a<l, 0<a<T, (51) 


with the initial condition 
Ug = u(0) = 0. (52) 


The exact solution for a = 1 is u(x) = 1+ V2 tanh(V2z + 5 log(3-4)); see 
[2]. 
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Table 4: Comparison between the approximation solution and exact solution 
of the presented method with TTM [2] for a = 1, T = 1, and N = 10 for 
Example 1 


x TTM [2] proposed method Exact 
0.0 0.0 0.0 0.0 
0.2 0.197773 0.197378 0.197374 
0.4 0.380422 0.379951 0.379949 
0.6 0.537449 0.537051 0.537050 
0.8 0.664285 0.664036 0.664037 
1.0 0.761671 0.761572 0.761594 


Table 5: Comparison between the absolute error of solution by our method 


with TTM [2] for a= 1 and T = 1, for Example 1 

x Error of proposed method Error of TTM [2] 
0.0 0.0 0.0 

0.2 1.4598 x 10-§ 7.2107 x 10-7 
0.4 1.5961 x 10-§ 1.7216 x 10-9 
0.6 4.6060 x 10-7 2.7186 x 10-3 
0.8 1.1006 x 10-® 3.3906 x 10-3 
1.0 1.9061 x 10-° 3.6117 x 1078 


In Figure 4, a comparison between the exact solution for a = 1 and the 
numerical solution for a = 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1 and T’ = 1 is 
shown. Also, Table 6 presents numerical solutions at some points of [0, 1] and 
for different values of a at T = 1. 


Table 7 shows the obtained maximum errors and ROC for a= 1, T = 1, 
and N = 5,10, 20,40, 80,160 Also, Figure 5 shows the numerical and exact 
solutions for a= 1, T =1, and N = 10. The numerical rate of convergence 
is highly consistent with our theoretical analysis results. 

Table 8 represents the present method and the achieved results of parti- 
cle swarm optimization (PSO) [2], modified homotopy perturbation method 
(MHPM) [2], Chebyshev wavelets (CW) [2], fractional variational itera- 
tion method (FVI) [2], Legendre wavelets method (LWM) [2], and Pad’e- 
variational iteration method (PVI) [2]. 

Table 9 presents a comparison between the exact solution for a = 1 and the 
numerical solution for T = 8. Also, Figure 6 shows a comparison between the 
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Table 6: Exact solutions and Numerical solutions of Example 2 for N = 10, 
T =1, and a = 0.8, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1 


a |0.3 0.6 0.7 0.8 0.9 0.95 0.99 0.999 1 Exact 

x 

0.1 |1.38 ]3.79 x 107! |2.65 x 107! ]1.92 x 107! |1.44 x 107+ [1.25 x 107! 1.13 x 107! ]1.10 x 107+ |1.10 x 107! }1.10 x 10+ 
0.2 |1.92 ]7.21 x 107! |5.25 x 107! |3.94 x 107! |3.04 x 107! |2.70 x 107-! |2.47 x 107! |2.42 x 107! |2.41 x 107! |2.41 x 107+ 
0.3 |2.14 1.02 7.80 x 107! |6.05 x 107! |4.82 x 107! |4.35 x 107! |4.02 x 107! 3.95 x 107! 3.95 x 107! ]3.95 x 1074 
0.4 12.24 1.29 1.02 8.22 x 107! 16.74 x 107! |6.16 x 10-! |5.76 x 10-! |5.68 x 10-1 |5.67 x 107! |5.67 x 10~! 
0.5 2.30 1.51 1.25 1.03 8.75 x 107! |8.10 x 107! |7.66 x 107! |7.56 x 107! |7.55 x 107! 17.55 x 107! 
0.6 12.34 1.70 1.45 1.24 1.07 1.01 9.64 x 107! |9.54 x 107! |9.53 x 107! |9.53 x 107! 
0.7 |2.36 1.84 1.62 1.43 1.27 1.20 1.16 1.15 1.15 1.15 
0.8 |2.37 1.96 LE 1.59 1.45 1.39 1.35, 1.34 1.34 1.34 
0.9 2.38 2.05 1.89 1.74 1.61 1.56 1.53 1.52 152) 1.52 
1.0 |2.39 2,12 1.99 1.87 1.76 1.72 1.69 1.69 1.68 1.68 

Table 7: Maximum absolute errors and ROC of Example 2 for a = 1, T = 1, 


and N = 5, 10, 20,40, 80, 160 


N Maximum Absolute Error ROC 
5 6.35 x 1078 - 

10 3.63 x 107° 7.45 
20 3.63 x 107° 3.32 
40 1.73 x 1077 4.39 
80 7.05 x 1079 4.62 
160 2.98 x 10719 4.57 


Table 8: Comparison of the numerical solutions of the equation in Example 
2 witha =landT=1 


z |SJOM [2] MHPM 2] [PSO 2] [CW P| |FVI P| [PVI 2] [LWM [2] [Our Method | Exact 

0.6 [1.007291 | 1.370240 [1.296320 [1.349150 [1.331462 [1.873658 |1.296302 | 0.953552 [0.953567 
0.7 |1.253674 | 1.367499 [1.416139 [1.481449 1.497600 [2.112944 [1.416311 | 1.152926 [1.152950 
0.8 [1.467499 | 1.794879 [1.506936 [1.599235 [1.630234 |2.260134 |1.506913 | 1.346363 |1.346365 
0.9 [1.629901 | 1.962239 1.569252 |1.705303 [1.724439 |2.339134 |1.569221 | 1.526897 _ |1.526913 
1.0 |1.787222 | 2.087384 1.605580 [1.801763 [1.776542 [2.379356 [1.605571 | 1.689487 [1.689500 
Table 9: Comparison between the exact solution and numerical solutions of 


Example 2 for a=1, T = 8, and N = 80 


Numerical solution 


Exact solution 


Error 


1.346362994 
2.246290755 
2.395782816 
2.412338083 
2.414131848 
2.414445422 
2.414746423 
2.415345681 
2.416609669 
2.418416749 


1.346363655 
2.246285959 
2.395756424 
2.412281528 
2.414012382 
2.414192625 
2.414211383 
2.414213335 
2.414213538 
2.414213559 


6.6128045 x 1077 
4.7957279 x 10~® 
2.6391922 x 107° 
5.6554231 x 107° 
1.1946588 x 1074 
2.527976 x 1074 
5.3504015 x 1074 
1.1323455 x 1073 
2.3961302 x 1073 
4.2031900 x 1073 
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2.5 1 1 


0.5 
*® = a=0.999 


<n a=1 
—*— Exact 
(0) 1 l l l l l l 1 
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Figure 4: Comparison between exact solution of Example 2 for a = 1 and 
numerical solutions for a = 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 1 and T = 1 


1.8 1 


—#—Numerical solution 
1.6;- ~ - - Exact solution 


0.4/—- 


nol 


Figure 5: Comparison between the exact solution and numerical solution of 
Example 2 for a=1, T=1, and N = 10 
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2.5 


—*- numerical solution 
- = =exact solution 


0 | ! | ! | | ! 


0 1 2 3 4 5 6 7 8 


Figure 6: Comparison between the exact solution and numerical solutions of 
Example 2 for a= 1, T = 8, and N = 80 


6 Conclusions 


This paper proposed a high-order compact finite difference method for the 
Riccati problem. The convergence analysis has been discussed. The numer- 
ical results presented in Tables 1-9 showed that the method is effective and 
that the numerical experiment is very consistent with our theoretical analysis 
results. 
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